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ABSTRACT 

We present the first fuhy 3D MHD simulation for magnetic channeling and con- 
finement of a radiatively driven, massive-star wind. The specific parameters are chosen 
to represent the prototypical slowly rotating magnetic O star 0^ Ori C, for which cen- 
trifugal and other dynamical effects of rotation are negligible. The computed global 
structure in latitude and radius resembles that found in previous 2D simulations, with 
unimpeded outflow along open field lines near the magnetic poles, and a complex equa- 
torial belt of inner wind trapping by closed loops near the stellar surface, giving way to 
outflow above the Alfven radius. In contrast to this previous 2D work, the 3D simula- 
tion described here now also shows how this complex structure fragments in azimuth, 
forming distinct clumps of closed loop infall within the Alfven radius, transitioning in 
the outer wind to radial spokes of enhanced density with characteristic azimuthal sep- 
aration of 15 — 20°. Applying these results in a 3D code for line radiative transfer, we 
show that emission from the associated 3D 'dynamical magnetosphere' matches well 
the observed Ha emission seen from 9^ Ori C, fitting both its dynamic spectrum over 
rotational phase, as well as the observed level of cycle to cycle stochastic variation. 
Comparison with previously developed 2D models for Balmer emission from a dynam- 
ical magnetosphere generally confirms that time-averaging over 2D snapshots can be a 
good proxy for the spatial averaging over 3D azimuthal wind structure. Nevertheless, 
fully 3D simulations will still be needed to model the emission from magnetospheres 
with non-dipole field components, such as suggested by asymmetric features seen in 
the Ha equivalent- width curve of Ori C. 

Key words: MHD — Stars: winds — Stars: magnetic fields — Stars: early-type — 
Stars: rotation — Stars: mass loss 



1 INTRODUCTION 

Massive, hot (OB type) stars lack the vigorous subsurface 
convection thought to drive the dynamo central to the mag- 
netic activity cycles of the sun and other cool stars. Nonethe- 
less, new generations of spectropolarimeters have directly 
revealed large-scale, organized (often significantly dipolar) 
magnetic fields ranging in strength from 0.1 to 10 kG in sev- 
eral dozen O B stars (e.g. [Donati et al.||2002l [20061 [Hubrig 
eFal. 2006; Petit et al. 2008; 'Grunhut et al.i2009| Martins 



et al.^^2010J with ongoing monitoring and surveys carried 
out by the MiMeS (Magnetism in Massive Stars) consor- 



tium dWade et aH2012p . 



A key theoretical issue, examined in a series of papers 
(lud-Doula & Owocki 2002; ud-Doula et al. 2008, 2009), is 



how such strong fields can, in combination with stellar rota- 
tion, channel and confine the radiatively driven stellar wind 
outflows of such OB stars, leading to the formation of a 
circumstellar magnetosphere. In relatively rapidly rotating 
magnetic Bp stars (with periods of order a day), magneti- 
cally trapped wind material above the Kepler co-rotation ra- 
dius is centrifugally supported, and so builds up into a stable 
centrifugal magnetosphere (iTownsend & Owocki 2005), of- 
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ten leading to rotationally modulated, high velocity Balmer 
line emission (Townsend et al. 2005 ). 

However, likely because of the angular momentum loss 
in their stronger magnetic winds, most 0-stars with detected 
surface magnetic fields are slow rotators^ with periods rang- 
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ing from 15 days for 0^ Ori C to 538 days for HP 191612 
and 55 years for HP 108 (|Stahl et al.||2008| [Howarth et"ar 



2007| [Martins et al. 2010), implying a very limited dy- 



namical effe ct of rotation on th e wind magnetic channel- 
ingM Indeed, Gagne et al. (2005) showed that non-rotating. 



two-dimensional (2P) magnetohydrodynamic (MHP) simu- 
lations match well the temperature, luminosity, and rota- 
tional phase-dependent occultation of the X-ray emitting 
plasma around 0^ Ori C. 

For the even more slowly rotating case of HP 191612, 



Sundqvist et al. (2012) have used similar 2P MHP simula- 



tions to show that both the strength and rotational phase 
variation of the observed Ha emission can be well repro- 
duced in terms of a dynamical magneto sphere] in this case, 
the higher 0-star mass loss rate means that, even without 
significant centrifugal support, the transient magnetic sus- 
pension and subsequent gravitational infall of magnetically 
trapped material yields a high enough circumstellar density 
to give strong Balmer line emission (see also Grunhut et al. 
2012 for a similar application to another 0-star HP 57682). 

To provide a sound basis for interpreting the Balmer 
emission in the growing sample of magnetic 0-stars being 



discovered in the MiMeS project ( |Wade et al.||2012| [Petit 
et al.||2012 ), the present paper applies this concept of a dy- 
namical magnetosphere to analyzing the rotationally mod- 
ulated Ha line emission of the prototypical magnetic 0-star 
6^ Ori C, but now using, for the first time, fully 3P MHP 
simulations of the magnetic channeling and confinement of 
its radiatively driven stellar wind. 

This extension to 3P represents a significant advance 
over the Sundqvist et al. (2012) approach, which was based 



on simulations that assume 2P axisymmetry, and so artifi- 
cially suppress any variations in the azimuthal coordinate 
(j). To mimic the expected azimuthal fragmentation of struc- 
ture in a 3P model, Sundqvist et al.|(|2012) computed spec- 
tra from multiple snapshots of a 2P MHP simulation with 
large stochastic variations, which were then time- averaged 
as a proxy for the radiation-transport spatial averaging over 
azimuth in a full 3P model. 

The 3P simulations here compute this lateral structure 
directly, applying the results to a full 3P radiative trans- 
fer model of the Balmer emission, and its variability with 
stochastic fluctuations of the wind, as well as with the chang- 
ing observer perspective over rotational phase. The next sec- 
tion (§2) describes the basic MHP method, along with the 
parameters, boundary conditions, and numerical grid used. 
This leads (§3) to a general description of the resulting 3P 
structure and evolution, with emphasis on the dense, near- 
equatorial, magnetically confined gas that makes up the dy- 
namical magnetosphere. These simulation results are then 
applied (§4) within full 3P radiative transfer calculations to 
derive the Balmer line emission and its variation. We con- 
clude (§5) with a summary discussion of the implications for 
understanding the circumstellar emission of magnetic mas- 
sive stars, along with an outline for future work. 



^ One exception is Plaskett's star ( [Grunhut et al.|2012| ), in which 
mass accretion from the close binary companion may have spun 
up the inferred magnetic star. 



2 3D SIMULATION METHOD 
2.1 Basic formalism 

For our 3P dynamical models we follow the basic methods 
and formalism developed for 2P simulations by |ud-Poul"a| 
& Owocki (2002), as generalized by Gagne et al. (2005) to 



include a detailed energy equation with optically thin radia- 
tive cooling (Mac Ponald Bailey||198T ). Since 3P model- 
ing, especially MHP, is computationally intensive, we now 
use the parallel version of the publicly available numerical 
code ZEUS-3P called ZEUS-MP (Hayes e t al][2006t , which 
is an independent code but employs essentially the same 
numerical algorithms, although in a parallel programming 
environment. 

As in |ud-Poula &: Owocki ( 2002[ ), the treatment of ra- 
diative driving by line-scattering follows the standard |Cas-| 
tor, Abbott & Klein (1975 hereafter CAK) formalism, cor- 
rected for the finite cone angle of the star, using a spher- 
ical expansion approximation for the local flow gradients 



( jPauldrach, Puis Kudritzki|1986| [Friend Abbott|1986 ) 
This ignores non-radial components of the line-force that 
may result from velocity gradients of the flow in the lateral 
directions. Although in purely hydrodynamic models these 
forces can have significant impact, e.g. by inhibiting wind 
compressed disk (Owocki et al. 1996|, in magnetic models 



they are generally small compared to the strong non-radial 



force components associated with the magnetic field ( ud- 
Poula||2003l . It also ignores the extensive small-scale wind 



structure and clumping that can arise from the intrinsic in- 
stability of line driving (| Owocki et al. 1988), since model- 
ing such dynamics requires a non-local line-force integration 
that is too computationally expensive to be tractable in the 
3P MHP simulations described here. 



2.2 Numerical grid 

The computational grid and boundary conditions are again 
similar to those used by ud-Poula & Owocki (2002), as gen- 
eralized by )ud-Poula et al. (2008) to include the non-zero 
azimuthal velocity at the lower boundary arising from rota- 
tion. Flow variables are specified on a fixed 3P spatial mesh 
in radius, co-latitude, and azimuth: The mesh 

in radius is defined from an initial zone i = 1, which has a 
left interface at ri = i?*, the star's surface radius, out to a 
maximum zone (i — rir — 300), which has a right interface 
at rsoo = lOi?*. Near the stellar base, where the flow gradi- 
ents are steepest, the radial grid has an initially fine spacing 
with An = 0.0048i?*, and then increases by 1% per zone 
out to a maximum of Ar299 = 0.095i?*. This is fine enough 
to resolve gradients in the subsonic base, but is larger than 
used in previous models; this allows for a larger time step, 
set as 0.3 of the Courant time, giving a typical time step of 
0.1 s. 

The mesh in co- latitude uses ne = 100 zones to span 
the two hemispheres from one pole, where the j — 1 zone 
has a left interface at Oi =0°, to the other pole, where the 
j = ne = 100 zone has a right interface at ^loo = 180°. 
To facilitate resolution of compressed flow structure near 
the magnetic equator at ^ = 90°, the zone spacing has 
a minimum of AO49 = A^so = 0.29° about the equator, 
and then increases by 5% per zone toward each pole, where 
AOi = A6>99 = 5.5°. 
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The mesh in azimuth is uniform, with = 120 zones 
of fixed size A(/) = 3° over the full 360° range. We also 
performed an initial test run with twice the number of az- 
imuthal zones for our future resolution study. But prelimi- 
nary results show that although there are some differences in 
finer structure during initial phases of breakup in azimuth, 
the asymptotic azimuthal structure scale is similar to that 
found here. 



2.3 Boundary conditions 

Boundary conditions are implemented by specification of 
variables in two phantom or ghost zones. At the outer ra- 
dius, the fiow is invariably super- Afvenic outward, and so 
outer boundary conditions for all variables (i.e. density, en- 
ergy, and the radial and latitudinal components of both 
the velocity and magnetic field) are set by simple extrap- 
olation assuming constant gradients. On the other hand, at 
the inner radius, where the surface of the star is defined, the 
fiow is sub- Afvenic. In general, the number of characteristics 
pointing into the grid determines the number of boundary 
conditions needed to be specified ( Bogovalov||1997 ) . For su- 
permagnetosonic fiow, this implies that all seven variables 
(density, internal energy or pressure, three components of ve- 
locity and two transverse components of the magnetic field) 
must be specified. But for sub-magnetosonic fiow here, only 
four characteristics point into the grid (entropy, fast, slow 
and Afvenic). Thus, four variables must be specified and 
the remaining three of the variables need to be extrapolated 
( Bogovalov|1997 ) . These conditions are briefiy described be- 
low. 

In the two radial zones below i = 1, we set the radial 
velocity Vr by constant-slope extrapolation, and fix the den- 
sity at a value po chosen to ensure subsonic base outfiow 
for the characteristic mass fiux of a ID, nonmagnetic CAK 
model, i.e. po ~ M/(47ri?^a/6). These conditions allow the 
mass fiux and the radial velocity to adjust to whatever is ap- 
propriate for the local overlying fiow (Owocki, Castor, and 
Rybicki 1988). In most zones, this corresponds to a subsonic 
wind outfiow, although infiow at up to the sound speed is 
also allowed. 

Magnetic fiux is introduced through the lower bound- 
ary as the radial component of a dipole field Br{R^,0) = 
Bocos{0), where the polar field Bo = 1100 G is adopted 
from the observationally inferred value for 0^ Ori C (Do- 
nati et al.| [2002). The latitudinal components of magnetic 
field. Be J and fiow speed ve, are again set by constant slope 
extrapolation. The azimuthal field 5^ is likewise set by ex- 
trapolation, while the azimuthal speed is set by the assumed 
stellar rotation, V(f){R^,0) — Vrot sin ^. 

While the code is fully 3D, the choice of spherical co- 
ordinates leads to numerical difficulties near the coordinate 
singularity along the polar axis, which is treated with re- 
fiecting boundary conditions. Because this prohibits any fiow 
across the axis, we ignore the slow, oblique rotation inferred 
for 0^ Ori C, since at a period of ^15 d, it is dynamically 
unimportant. However, to avoid pockets of rarefied gas near 
the stellar surface that lead to very small Courant time, we 
do introduce a small dynamically insignificant field- aligned 
rotation with an equatorial surface speed Vrot = lOkms"^. 



Table 1. Stellar and wind parameters 



Name 


Parameter 


Value 


Mass 




4OM0 


Radius 




8Rq 


Luminosity 


L* 


1.5 X IG^Lo 


CAK a 


a 


0.5 


Line strength 


Q 


700 


Mass loss rate 


M 


3.3 X lO-'^Moyr-i 


Terminal speed 




3200 km s-i 


Polar field 


Bpole 


1100 G 



2.4 Stellar parameters and initial condition 

The adopted stellar and wind parameters are summarized in 
Table 1. Using these parameter values, we first relax a non- 
magnetic, spherically symmetric wind model to an asymp- 
totic steady state; this yields the quoted mass loss rate and 
terminal speed, which are in very good agreement with the 
standard CAK scaling, adjusted for a finite-disk and non- 
zero sound speed (Owocki & ud-Doula 2004). This relaxed 
CAK steady state is used for the initial condition in density 
and velocity. The temperature is initially set to the associ- 
ated stellar effective temperature, Teff ^ 40, 000 K, which 
also subsequently sets a fioor for the temperature, to ap- 
proximate a radiative equilibrium in which photoionization 
heating counters the cooling by line emission and adiabatic 
expansion. 

The magnetic field is initialized to have a simple 
dipole form with components Br = B o{Ri>. / r)^ cos ^ Be — 
{Bo/2){R^/rfsmO, and = 0, with Bo the polar field 
strength at the stellar surface. From this initial condition, 
the numerical model is then evolved forward in time to study 
the dynamical competition between the field and fiow, as de- 
tailed in ^3. 



Following ud-Doula & Owocki (2002), the stellar, wind, 
and magnetic parameters in Table 1 can be use(|^to define 
a dimensionless "wind magnetic confinement parameter" , 



MVoo 



14, 



(1) 



where for a dipole, the equatorial field is just half the polar 
value, Beq = Bo/2. Since 77* > 1, the magnetic field domi- 
nates the wind outfiow near the stellar surface. But because 
the magnetic energy density has a much steeper radial de- 



cline {B^ 



for a dipole) than the wind {pv 



if near terminal speed), the wind can strip open field lines 
above a a characteristic Alfven radius, set approximately by 
dud-Doula et al.|2008t 

- 0.3 + 77^^^2.23. (2) 



R. 



^ We emphasize here that mass loss rate and flow speed used 
to estimate this magnetic conflnement parameter (and thus the 
associated Alfven radius) should not, as is sometimes done, be 
based on any observationally inferred values using non-magnetic 
analysis for the star in question, as the magnetic fleld may signif- 
icantly influence the predicted circumstellar density and velocity 
structure. Rather, 77* is deflned using the mass loss and speed ex- 
pected for a corresponding non-magnetic star of that spectral and 
luminosity class, as inferred from scaling laws (e.g., |v'ink et aE] 
2001 1 based on line-driven-wind theory. 
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Figure 1. A near stellar surface view of an iso-density (logp = 
— 12.5 in cgs units) surface, colored with radial velocity Vr, ranging 
from — lOOkms"-*^ (red) to +100 km s"-*^ (blue). Note how some 
of the gas material within the closed loop region falls back onto 
the stellar surface, while material in the open field (polar) region 
streams radially outward. The slowly moving (or infalling) region 
depicted as white (red) color indicates the lateral extent of the 
magnet osphere. The circle (cyan) represents the stellar surface. 



3 MHD SIMULATION RESULTS 
3.1 Snapshot of 3D structure 

As background for the Balmer emission line computations 
in §4, let us here characterize the spatial structure and time 
evolution resulting from our 3D MHD simulation. 

The characteristic time for the CAK wind here to flow 
from the sonic point near the stellar surface to the outer 
boundary at lOi?* is about tfiow ~ 25 ks. For MHD sim- 
ulations we find in practice (ud-Doula & Owocki |2002 ) it 
takes ~ 10 — 15 such flow times for resulting structure to 
reach a state that is independent of the initial conditions. 
To provide an extended sampling of random variations in 
this asymptotic state, we run the 3D simulations here to a 
final time t = 1000 ks, corresponding to roughly 40 wind 
flow times. 

Figure [l] shows a close-in, equatorial view of the di- 
chotomy between open field over the poles and closed loops 
around the equator, with density rendered as an iso-density 
surface for logp = —12.5 g/cm^. The color represents the 
associated radial (in grid coordinate system) component of 
the flow velocity, ranging from -100 kms~^ inflow (red) to 
+ 100 kms""*^ outflow (blue). Note that the polar regions 
have mainly just outflow, reflecting the open nature of the 
field. In contrast, the equatorial regions have a complex com- 
bination of outflow and inflow, reflecting the trapping of the 
underlying base wind within closed loops, with subsequent 
gravitational fall-back onto the stellar surface. 

A similar pattern of wind trapping and subsequent infall 
was also seen in 2D simulations (ud-Doula & Owocki 2002] 
ud-Doula et al.|2()08 ), manifested there by a complex 'snake' 
pattern extending over radius and latitude, but formally co- 
herent in azimuth within the imposed 2D axisymmetry. This 
3D simulation now shows how, because of spontaneous sym- 



metry breaking, this infall becomes azimuthally fragmented 
into multiple, dense, infalling clumps, distributed in both 
latitude and azimuth on a scale that is comparable to the 
latitudinal scale of the 2D snakes. 

To provide a more complete picture of this complex, 
3D wind structure, figure [5] shows both equatorial (top) and 
polar (bottom) views for two values of the iso-density sur- 
face, namely log p = —12.5 (left) and log p = —13 (right), 
again with the same color coding for radial velocity. The po- 
lar view now shows clearly the azimuthal organization and 
scale of the trapped equatorial structure, and how it be- 
comes extended outward into radial "spokes". While com- 
plex, the lower iso-density plot on the right shows a general 
trend from dominance of infall in the inner region (marked 
by red) to outflow at larger radii (marked by blue). 

Overall, this represents a 3D generalization of the dy- 
namical magnetosphere concept developed to characterize 
2D simulations, and their application to modeling Balmer 



line emission (Sundqvist et al. 2012). The essential fea- 



ture is that the trapping and compression of wind mate- 
rial in closed loops results in a substantial enhancement of 
the time-averaged density near the magnetospheric equator, 
leading then to an associated enhancement in line emission, 
as we demonstrate quantitatively in §4. 



3.2 Time development and evolution of wind 
structure 

Let us next consider how such complex spatial structure de- 
velops from the steady, spherical initial condition, and how 
it subsequently evolves into a stochastic variation associ- 
ated with wind trapping, breakout, and infall. Since much 
of the key structure is near the magnetic equator, let us fol- 
low the approach developed by ud-D oula et al.| ( |2008| ) that 
defines a radial mass distribution drue/dr within a speci- 
fied co-latitude range (taken here to be 10°) about the 
equator. 



dme{r, 0, t) 
dr 



: 27rr' 



I 



7r/2 + A6'/2 



/2-A6'/2 



p(r, ^, 0, t) sin^ 



(3) 



In previous 2D models, the imposed axisymmetry means this 
dnie/dr is a function of just radius and time, and so can be 
conveniently rendered as a colorplot to show the time and 
spatial variation of the equatorial material. In the present 
3D case, there is an added dependence on the azimuthal 
coordinate 0, but upon averaging over this, we again obtain 
a function of just radius and time. 



dme{r, t) 
dr 



27r 



drrie 
dr 



(r, 0, t) dcj) . 



(4) 



Figure |3] compares the time and radius variation of the 
dfhe/dr from the 3D simulation (left) with the correspond- 
ing drrie/dr from an analogous 2D model (right) with iden- 
tical parameters. The dashed horizontal line, marking the 
estimated Alfven radius Ra ~ 2.23i?*, roughly separates 
the regions of outflo^^0 along open field lines above Ra^ 

^ Actually, for the 2D case, the separation is more distinct in 
isothermal models shown by |ud-Doula et al.| ( |2008 ^; in the 2D 
simulations here, the inclusion of a full energy balance leads to 
episodes of field reconnection that pull relatively high-lying ma- 



isodensity 1 0""^ ^ g/cm^ ^ 
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Figure 2. Iso-density surfaces at log(p)=-12.5 (left) and log(p)=-13.0 (right) for a fixed time snaphot t = 1 Ms, viewed from the equator 
(top panels) or pole-on (bottom panels). As in figure [l] the color represents radial velocity. The equatorial 'belt' visible as white band 
(low velocity) represents the latitudinal extent of the magnet osphere. Note that substantial material within that region is falling back 
onto the stellar surface, thus forming the dynamical magnetosphere. The cyan circle again represents the stellar surface. 



from the complex pattern of upflow and downflow for mate- 
rial trapped in closed magnetic loops below Ra- In the 2D 
models, the complex infall pattern occurs in distinct, snake- 
like structures, but in the 3D models, these break up into 
multiple azimuthal fragments that, upon averaging, give the 
3D panel a smoother, more uniform variation. More impor- 
tantly, note the clear enhancement in the mass just below 
the Alfven radius, representing what we characterize here as 
the dynamical magnetosphere. 

To illustrate the time evolution of the azimuthal 
breakup and variation in 3D wind structure, the upper pan- 
els of figure |4] use a "clock format" to show the time progres- 
sion of an arbitrary 36° segment of dme/dr{r^(j)^t) at time 
steps of 50 ks, ranging over t = 50 — 500 ks for the left panel, 
and continuing through t = 550 — 1000 ks in the right panel. 
Following the spherically symmetric initial condition, the 
magnetic field first channels the wind into a compressed, az- 
imuthally symmetric, disk- like structure (t = 50 & 100 ks). 



terial back to the surface. The reasons for these subtle differences 
will be investigated in a follow-up study. 



which however then breaks up (by t = 150 ks) to a series 
of radial spokes with an azimuthal separation of roughly 
A(/) 15 — 20°. Remarkably, following this initial breakup, 
the azimuthal structure keeps the same overall character and 
scale within a complex stochastic variatiorj^ 

Note in particular that the structure over the full 360 ° , 
as shown for the final time t — 1000 ks in the lower row 
of figure [4] is quite similar to that shown for the evolving 
segments in the right-side upper row for t = 550 — 1000 ks. 
Moreover, the ^ 500 ks lower bound for this settling time 
corresponds to the time in figure [3] for the formation of the 
enhanced dirie/dr that characterizes the dynamical magne- 
tosphere. 

As noted in ^3.1| the inferred azimuthal scale of 15 — 20° 



Initial test runs done at twice the azimuthal grid resolution 
from that used here show an initially smaller scale for the breakup, 
which however subsequently (by ~ 300 ks) settles into a scale that 
is similar to what is found here. Details will be given in a future 
paper that carries out a systematic resolution study for such 3D 
MHD simulations. 
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Figure 3. Equatorial mass distribution drue/dr (in units 
Mq/R^) for the azimuthally averaged 3D model (left) and for 
a corresponding 2D model (right), plotted versus radius and 
time, with the dashed horizontal line showing the Alfven radius 
Rx ~ 2.23i?*. Following a similar adjustment to the initial con- 
dition, the long-term evolution of the 2D models is characterized 
by a complex pattern of repeated strings of infall, whereas the 
3D azimuthally averaged model settles into a relatively smooth 
asymptotic state characterized by an enhanced mass near and 
below the Alfven radius. 



is comparable to the scale of infalling "snakes" found in 2D 
simulations, which typically extend over some fraction of the 
loop length between closed footpoints. For the present case, 
the maximum latitude of a closed loop is about 45°, giving 
a typical angle extent of about 15 — 20° for both the snakes 
and the azimuthal breakup. Future studies should examine 
whether these structure scales do in fact follow the expected 
loop closure latitude, which should scale with magnetic con- 
finement parameter as arcsin?77"'^^^. 

In summary, the equatorial structures seen in previous 
2D models now become azimuthally fragmented with a char- 
acteristic separation of A(/) 15 — 20°, which however upon 
averaging in azimuth, have a radial and latitudinal variation 
that is similar to what is obtained from time- averaging in 
2D models. This has important implications for modeling 
the associated Balmer line emission, as we discuss next. 




Figure 4. Time evolution of dnie/dr for an arbitrary 36 ° az- 
imuthal segment of the 3D model, at marked time intervals of 
50 ks ranging over f = 50 — 500 ks in the upper left panel, extend- 
ing to t = 550 — 1000 ks in the upper right. The bottom panel 
shows the full variation over 360 ° for the final simulation time 
t = l Ms. As in figure [3] the units of dme/dr are Mq/R^. 



a lower boundary condition, and solve the equations of radia- 
tive transfer only in the circumstellar structure. Hydrogen 
NLTE departure coefficients are estimated from a spheri- 
cally symmetric unified model atmosphere calculation (us- 
ing FASTWlND, |Puls et al.|2005] ). Moreover, since the energy 
equation in the MHD simulation gives only a rough esti- 
mate of the wind temperature structure, we also take this 
from the FASTWIND model, except for in shock-heated re- 
gions with T > 10^, where we set the Ha source function 
and opacities to zero. For further discussion about these and 
other assumptions entering the Balmer line profile calcula- 



tions, see Sundqvist et al. (2012) and references therein 



4 BALMER LINE EMISSION: MODEL VS. 
OBSERVATIONS 

Now let us examine to what extent the 3D MHD wind simu- 
lations presented here can reproduce the observed Ha emis- 
sion and variability in Ori C. 

Applying the 3D radiative transfer method developed 



by Sundqvist et al. (2012| to the 3D MHD simulations de- 
scribed above, we compute synthetic Ha profiles by solving 
the formal integral of radiative transfer in a cylindrical co- 
ordinate system aligned toward the observer. The angle a 
defines the observer's viewing angle with respect to the mag- 
netic pole, and, for given inclination i and obliquity /3, 



cos a = sin i sin /3 cos $ + cos i cos /3 



(5) 



gives the observer's viewing angle as a function of rotation 
phase Due to the slow rotation of the targeted star, we 
may neglect any dynamical effects of rotation and use the 
same simulation for all obliquity angles /3. The computa- 
tions further use a pre-specified photospheric line profile as 



4.1 Periodic rotational phase variability 

Figure [5] displays synthetic Ha line-center surface bright- 
ness maps and flux profiles viewed from above the mag- 
netic pole and equator, corresponding approximately to ro- 
tational phases of 0.0 and 0.5, respectively, according to 
the ephemeris of |Stahl et al.| ( [2b08| . As in|Sundqvist et al. 
( |2012| ), we interpret the Ha variability within the frame- 
work of a dynamical magnet osphere, but now using the full 
3D dynamical model for the azimuthal wind structure. For 
the pole-on view, the surface brightness of line-center emis- 
sion in the lower panel of figure [5] follows closely the lateral 
fragmentation seen from the corresponding iso-density plots 
in the bottom row of figure [2] Similarly, there is a good cor- 
respondence for the equatorial view, but now the projection 
leads to an overall lower level of the total emission, associ- 
ated with the smaller projected surface area of the struc- 
ture. As shown by the right panels of figure [5] this leads 
to a stronger level of emission in the line profiles for the 
polar view. But note that even for the equatorial view, the 
circumstellar emission still partially refills the photospheric 
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Figure 5. Synthetic Ho; surface brightness maps (left panels) 
and flux profiles (right panels) for observer viewing angles above 
the equator (upper panels) and magnetic pole (lower panels), as 
calculated from a snapshot of the 3D MHD wind simulation taken 
at 1000 ksec after initialization. 
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Figure 6. Observed and synthetic Ha dynamic spectra of 0^ Ori 
C. Synthetic spectra calculated from patching together 2D model 
snapshots (panel 1), a single 3D model snapshot (panel 2), and 
from a random selection of 3D snapshots (panel 3). The vertical 
small lines in panel 4 represent the velocity range of the nebular 



emission. 



absorption profile. Most of the emission arises from the rel- 
atively dense material trapped in the closed loops near the 
stellar surface, with only about 10-20% arising from the ex- 
tended radial spokes. Overall, the magnetospheric structures 
found here have quite different morphology from compres- 
sive structures generated by instabilities of line-driving in 



non-magnetic stars ( Feldmeier et al. 



2003). 



1997 



Dessart & Owocki 



Within the context of the MiMeS project, 18 high- 
resolution, high-SNR observations of 0^ Ori C were obtained 
with the spectropolarimeter ESPaDOnS at the Canada- 
France-Hawaii telescope (4 from [Petit et al.| ([2008| and 14 



from MiMeS) which were used to produce the dynamical 
spectra of Ha shown in figure [6] Leveraging the high spec- 
tral resolution (4 kms~^), we remove the strong nebular 
line by fitting a polynomial curve to the adjacent sides of 
the line profile. This procedure provides us with an indica- 
tion of how smoothly these lines vary in phase. Nonetheless, 
localised features within this velocity range (illustrated by 
small vertical lines in the fourth panel in figure [gJ could still 
be hidden. 

Figure [g] compares observed dynamic Ha spectra of 0^ 
Ori C over its rotational phase (far right) with 3 models that 
each assume i = /3 = 45°, but which have increasing levels 
of sophistication progressing from left to right. 

The first panel is computed from 100 randomly selected 
2D simulation snapshots that have been patched together to 
form an 'orange slice' pseudo-3D model with an azimuthal 
coherence of 3.6°. This patch approach is similar to that 
used to reproduce periodic Ha variability in the magnetic 
0-stars HD 191612 and HD 57682 from 2D M HD simulations 
dSundqvist et al.|2012| [Grunhut et al.|2012 ). 

The second panel then displays spectra computed from 
a single snapshot of the 3D simulation, and illustrates that 
such full 3D model spectra are in fact qualitatively very 
similar to the averaged 2D models, as was suggested by 



Sundqvist et al. (2012) 



Finally, the third panel in Figure [6] is constructed from 
a random selection of 3D snapshots with phases that co- 
incide with the observations seen in the fourth panel; it is 
discussed further in the following subsection. 

Overall, all three models give reasonable and compara- 
ble agreement with the observed dynamic spectrum. This 
suggests that one can indeed use the more simplified ap- 
proaches to explore how observational diagnostics constrain 
the properties of a dynamical magnet osphere, calibrated 
with selected comparisons with more computationally de- 
manding, and physically realistic 3D models. 

Although the simulations reproduce quite well both the 
magnitude and phase of the observed rotational variation, 
the observed profiles exhibit an asymmetry about phase 0.5 
that is not found in our numerical models; such asymmetries 
may have their origin in, e.g., non-dipolar components of 
the magnetic field, consideration of which we defer to future 
study. Moreover, to match the absolute level of emission, we 
have increased the Ha scaling invariant y^fciM / (vooR*)^^'^ 
of the underlying simulations by a factor of 2.5, where /d 
is an overdensity correction factor that accounts for small- 
scale clump structure (e.g., Puis et ar]|2008 ). The need for 
this increase may indicate the MHD simulations either as- 
sume a somewhat too low mass- loss rate for 0^ Ori C, or 
underestimate the overdensities of the emitting Ha material 
around the magnetic equator. 



4.2 Stochastic short time-scale variabiHty 

The Ha emission of Ori C has been monitored regularly 
for nearly two decades and an extensive set of equivalent 
width measurements has been compiled by Sta hl et al.| 
( [2008 ) shown as the black dots in figure [t] phased accord- 
ing to the 15.4 d rotation period using the ephemeris of 
Stahl et al. (2008). Although the level of both emission 



and its modulation with rotational phase are very stable, 
there is some significant random variability present during 
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Figure 7. Observed and synthetic Ho; equivalent- width curves. 
The synthetic measurements are calculated from 100 random 3D 
simulation snapshots phased randomly. The quoted cr's measure 
standard deviations from functional fits (dotted green lines) to 
simulations and data. The blue squares denote the additional ob- 
servational data points that correspond to the dynamical spectra 
in figure [6] See text. 



all phases. Such stochastic Ha scatter may originate from 
small differences in the amount of material trapped in the 
dynamical magnetosphere at any given time. 

To examine the level of such short time-scale variabil- 
ity in our 3D MHD model, we randomly select simulation 
snapshots between 800-1000 ksec and from these calculate 
line profiles for random rotational phases. The results show 
that the model closely reproduces both the systematic vari- 
ation with rotational phase, as well as the stochastic devi- 
ation in the mult i- year data from the mean at each phase. 
The quoted standard deviations (Jsim ~ cTobs ~ 0.25 A have 
been estimated by calculating the variance from separate 
functional fits f{a) = a + 6|cosa| to the model and ob- 
servational data. (Due to the asymmetry, data are fit only 
between phases and 0.5.) Though the interpretation of 
this simple estimate is complicated by the fact that non- 
magnetic 0-stars with significant wind strength also show 
stochastic Ha variability (Markova et al.||2005), it never- 



theless suggests that in Ori C the observed Ha scatter 
is likely dominated by processes associated with the star's 
dynamical magnetosphere. 

This approach of randomizing variations in model 
equivalent width can also be used to account for the ef- 
fect of stochastic wind variations on the synthetic dynamic 
spectrum, as illustrated in the third panel of figure [6] We 
have aligned synthetic model phases to those for the subset 
of times with available spectra from our MiMeS observa- 
tions (blue squares in figurej7|. Comparison of the third 
and fourth panels of figure [Gjagain show an overall good 
agreement between the models and observations. 



5 DISCUSSION, CONCLUSIONS AND 
FUTURE WORK 

A central advance of the work described here is the incorpo- 
ration of all three spatial dimensions in the MHD simulation 



model. This allows us to follow directly the spontaneous 
breaking of the azimuthal symmetry that was imposed in 
previous 2D models, while still confirming a global struc- 
ture that is quite similar to that found in those models. At 
high latitude this is characterized by fast wind outflow along 
open field lines. At low latitudes, an equatorial belt of closed 
loops traps the wind into a complex pattern of upfiow and 
infall near the star, but this gives way again to open field 
outflow above the Alfven radius. 

The addition of a third dimension along the azimuth 
means that this complex pattern now fragments into dis- 
tinct regions of upflowing wind and infalling clumps. This 
leads to an associated fragmentation of the wind breakout 
above the Alfven radius, resulting in an alternating pattern 
of radial spokes of slow, high-density outflow, with a char- 
acteristic azimuthal separation of 15 — 20°, divided roughly 
equally between the dense spokes and a more rarefied flow 
between them. Such a separation angle is well above the 
adopted azimuthal grid size A(/) ^ 3°, and indeed a higher 
resolution test run with half this grid size gives asymptotic 
wind structure at a similar scale. As such, this scale does 
not appear to be strongly dependent on the grid resolution, 
though further investigation will be needed to test this and 
clarify what physical processes set the structure size. 

An additional focus here is the application of these 3D 
MHD simulation results to modeling the Ha line emission 
from the prototypical magnetic massive star 0^ Ori C, rep- 
resenting a 3D extension of the dynamical magnetosphere 
paradigm recently introduced by |Sundqvist et al.| ( |2012| ). As 
a proxy for the spatial averaging over the azimuthal vari- 
ations expected from a full 3D model, this previous work 
computed time-averaged spectra from a strictly 2D axisym- 
metric simulations with stochastically varying wind struc- 
ture. As demonstrated by the close similarity of the 2D vs. 
3D models for the dynamic spectra in figure [6] a key result 
here is to confirm a general validity for this proxy approach. 
This suggests then that future modeling of such rotational 
variation of line emission in the growing number of slowly 
rotating magnetic 0-type stars could be achieved with rea- 
sonable accuracy using much less computationally expensive 
2D MHD simulations to characterize the dynamical mag- 
netosphere. Of course, such an approach is only possible 
for the simple, axisymmetric field configuration of a mag- 
netic dipole. To account for deviations from a dipole, such 
as might be inferred from the asymmetric phase variation in 
the equivalent width light curve here, one will again need to 
carry out fully 3D simulations. 

Apart from this modest non-dipole asymmetry, both the 
2D and 3D dipole-based models here match quite well the 
rotational phase variation of Ha. But a further quite remark- 
able result of the present study is how well the synthetic 3D, 
"random time" model (third panel in figure [g]) reproduces 
also the observed stochastic variation in the Ha equivalent 
width, as demonstrated in figure [71 This suggests that the 
number and sizes of distinct structures generated in the sim- 
ulated dynamical magnetosphere may be quite comparable 
to that in the actual star. 

On the other hand, the magnetic 0-star HD 191612, 
which has a stronger magnetic confinement (77* ^ 50), shows 
a distinctly lower level of such random variation in Ha equiv- 



alent width (Howarth et al. 2007). This suggests that the 



general number and size of structures may depend on, e.g.. 
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the size or rigidity of the magnet osphere. This again requires 
further investigation through a systematic parameter study. 

Finally, while the overall agreement between 3D model 
predictions and observed Ho; variations is quite good, this 
is achieved with a modest adjustment in the assumed mass 
loss rate and/or wind clumping factor. Future work should 
examine the role of the line-deshadowing instability in in- 
ducing clumping in the wind acceleration region, and how 
well the numerical grid resolves this and other small-scale 
flow structures arising from wind compression in closed mag- 
netic loops. In addition, future 3D models should consider 
both dipolar and non-dipolar field topologies that are not 
aligned to the rotation axis. In summary, while the results 
here have provided an encouraging first step, there remains 
much work to develop realistic 3D models of the magneto- 
spheres and channeled wind outfiows from magnetic massive 
stars. 
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